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Using a recently developed method for calculating se- 
ries expansions of the excitation spectra of quantum lattice 
models, we obtain the spin-wave spectra for square lattice, 
5 = 1/2 Heisenberg-Ising antiferromagnets. The calculated 
spin-wave spectrum for the Heisenberg model is close to but 
noticeably different from a uniformly renormalized classical 
(large-S) spectrum with the renormalization for the spin- wave 
velocity of approximately 1.18. The relative weights of the 
single-magnon and multi-magnon contributions to neutron 
scattering spectra are obtained for wavevectors throughout 
the Brillouin zone. 



Thermodynamic properties and excitation spectra 
of two-dimensional quantum antiferromagnets have at- 
tracted much attention, cspcciaUy because of their po- 
tential relevance to high temperature superconductivity 
in the cuprate perovskitesB Methods based on high-order 
series expansions about the Ising model have proven to be 
very successful in accurate calculations of the thermody^ 
namic properties of two-dimensional antiferromagnetso'u 
One major limitation of these methods has been their in- 
ability to deal directly with dynamical properties or ex- 
citation spectra. These quantities have so far been stud-, 
ied, within the series approach, via, frequency momentsu 
and single mode approximations The reliability of us- 
ing only a few moments to reproduce spectral lineshapes 
has come into question, for example, in the case of the 
two-magnon Raman spectra.cl Currently efforts are un- 
derway to use a large number of frequency moments to 
obtain the dviiamical properties by numerical analytical 
continuation.Ll |-| 

Recently one of usa has shown that a series expansion 
method can be used to directly calculate excited state 
properties of quantum many-body systems. Here, we ap- 
ply this method to calculate the energies and spectral 
weights of the elementary excitations in square lattice, 
S = \/2 Heisenberg-Ising magnets by expansion around 
the Ising limit to order (J_l/Jz)^°. 

Let us briefly summarize the results. The spin-wave 
spectrum has an energy-gap for Jj^ < J^, which vanishes 
as the Heisenberg limit (J^ = Jz) is reached. Using se- 
ries extrapolation methods to estimate the dispersion in 
the Heisenberg limit, we find a spin- wave velocity which 



agrees with 1/5 expansions and other previous estimates. 
Along the line — Qy, the dispersion is nearly uniformly 
renormalized with respect to its classical value. Measur- 
able deviations are found in other directions. In partic- 
ular, along the line q-x + Qy = the spin- wave energy is 
maximized at (7r/2, tt/2) and it exceeds its value at (tf, 0) 
by about 7%. 

Another quantity that has been calculated by this 
method is the single-magnon contribution to neutron 
scattering, that is, the coefhcient associated with the 
delta-function in the dynamic structure factor. When 
compared with the equal-time transverse correlation 
function, it yields the relative weight of multi-magnon 
and single-magnon contributions to the neutron scatter- 
ing spectra at different wavevectors. We find that over 
substantial parts of the Brillouin zone approximately 20% 
of the total spectral weight is associated with the multi- 
magnon excitations even though the relative weight of 
such excitations vanishes near (0,0) and (7r,7r). 

The Heisenberg-Ising Hamiltonian under consideration 
is defined by 



W = J. 



(1) 



where the sum runs over nearest-neighbor pairs on a 
square lattice for which the lattice constant is the unit 
length (a = 1), and a = Jj_/ Jz- In the Ising limit, a — 0, 
there are two degenerate ground states and the single- 
magnon excitations are single spin-flips with respect to 
the Neel states. In this limit the excitations are purely 
local, or, alternatively, one could say that the magnon 
energies are degenerate over the entire band. For a ^ 0, 
the single-magnon states evolve into a set of states with 
a nonzero dispersion. The key to calculating the spin- 
wave dispersion is to construct an effective Hamiltonian 
for the states which are the natural, perturbatively con- 
structed extensions of the single spin-flip states at finite 
a. The effective Hamiltonian is then readily diagonalized 
by Fourier transformation. 

Because the spins on the two sublattices are oriented 
in opposite directions, the single spin-flip states naturally 
divide into two sets, those corresponding to S*^ = -1-1 and 
those with S*^ — —1. The effective Hamiltonian, which 
conserves S^, thus connects only the basis states with 
spin-flips on the same sublattice: this ensures that the 
spin-wave spectrum is degenerate between wavevectors 
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ilx, Qy) and {tt—Qx, n—qy). The full effective Hamiltonian 
to order a^'^ is presented in Table |[ 

For a 7^ 1 there is a gap in the spectrum, the minimum 
being at (0,0) and (7r,7r). The expansion for the gap is 

2 - (5/3)a2 + 0.317129630"* - 0.41923376a^ 
+ 0.27099699a* - 0.389433510*° + • • • (2) 

which agrees completely with the "mass gap" calculated 
by Zheng et alJ3 General arguments tell us that the gap 
must close in the Heisenberg limit, a = 1. Moreover, we 
expect that for small <? = |q| the spectrum has the form 
e(q) « [A{a) + B{a)q^]^/^, where A{a) ^ as a ^ 1. 
The spin-wave velocity in the Heisenberg limit is given by 
B{iy/^. In order to calculate the spin-wave velocity, we 
expand e(q) in powers of q, e{q) — C{a) + D{a)q^ + ■ • • , 
and identify C = A^/"^ and D = B/2A^/'^. Thus the 
square of the spin- wave velocity for the Heisenberg model 
is given by the a — > 1 limit of the series 

2C{a)D{a) = 4a2_2.305555a'* + 2.410512a*' 

- 3.064895a* + 4.100549a*" + ---. (3) 

Since we expect only a weak, energy-like singularity in 
this series at a = 1, we can sum it by Fade approxi- 
mants. The near-diagonal approximants [2/2], [2/3] and 
[3/2] give estimates of 2.779, 2.785 and 2.785 respectively, 
which leads to values for the spin- wave velocity c/ Ja of 
1.667, 1.669 and 1.669. In the largc-S* limit c/Ja = V2, 
so the quantum renormalization of the spin-wave veloc- 
ity is Zc ~ 1.18. This number is in excelleiit agreement 
with the high-order pcin-wave calculation^ and previ- 
ous indirect estimatesE2l using the hydrodynamic relation 
— Ps/Xi^i where ps is the spin-stiffness and x± the 
uniform transverse susceptibility. 

Away from the gapless points (0,0) and {tt,tt) the 
spin-wave spectrum for the Heisenberg model can be es- 
timated by direct Fade approximants. Along the line 
1x = lyi the spectrum, within our numerical uncertain- 
ties, is uniformly renormalized with respect to the clas- 
sical (large-S") spectrum. However, along the line qx = 
appreciable differences appear. In particular, the excita- 
tion energy at (tt, 0) appears to be a shallow local mini- 
mum along this line, and is lower than that at {tt/2, 7r/2) 
by about 7%. Note that at the classical and 1/5* lev- 
els the spin-wave spectrum is degenerate at (tt, 0) and 
(7r/2,7r/2). A plot of the dispersion relation along se- 
lected directions is shown in Fig. ^ A hint of the de- 
viations from classical spin-wave theory which we have 
found near (vr, 0) can be seen in the results of the projepj 
tor quantum Monte Carlo calculations of Chen et al.,i3 
even though those authors emphasize the similarity of the 
5 = 1/2 spin- wave spectrum and the uniformly renormal- 
ized classical spectrum. The first calculation of the spec- 
trum by spifcwave expansion to order 1/5*^ by Igarashi 
and Watabal3 also suggested that (tt, 0) would be a local 
minimum for the spin-wave spectrum; however, that min- 
imum is no longer presetit in the more recent spin-wave 
calculations by Igarashi. B 



We now turn to the spin-wave spectral weights. The 
magnetic neutron scattering cross section is proportional 
to the dynamic structure factor, given by the expression 

S'(q,cj)= f dte-"^*^e''i-''(5^(O,0)S'^(r,t) 

+sy{o,o)sy{r,t)). (4) 

We consider the T — limit, where the angular brackets 
refer to ground state expectation values. In general, we 
expect 5'(q, w) to consist of a sum of two parts, 

Siq,Lu)^Aiq)SiLO~eiq))+B{q,Lu). (5) 

Here e(q) is the spin-wave dispersion and A{q) is the 
residue or spectral weight associated with the spin- waves. 
B{q,uj) is associated with multi-magnon excitations, 
which are present because a single spin-flip in an anti- 
ferromagnet cannot be exactly represented as a superpo- 
sition of single-magnon states. Integrating S{q,u!) over 
all frequencies yields the equal-time correlation function 

S'(q) = J2 e''''''(5^(0, 0)5*=^ (r, 0) + 5^(0, 0)5*^ (r, 0)). 

r 

(6) 

To determine the residue A(q), we need to restrict the 
intermediate states that arise in the calculation of the ex- 
pectation values to be single magnon states, which gives 
the expression 

A{q) = J2 e"' ""(^"(0, 0)7'5"(r, 0) + 5^(0, 0)7^5^ (r, 0)), 

r 

(7) 

where V is the projection onto the manifold of single- 
magnon states. In the cluster expansions, each single- 
magnon state evolves with the coupling a, and is of the 
form 

1^.) = i^>+E^^»(")i")' (8) 

n 

where |i) is a single-magnon state in the Ising limit and 
\n) represents basis states (eigenstates in the Ising limit) 
which are not degenerate with the single-magnon states. 
However, the states {ipi) for different i are not orthogonal 
to each other when a 7^ 0. Thus in order to construct the 
projection operator, we need to define the overlap matrix 
gij =< ipi\ipj >■ Then the projection operator onlfl the 
single-magnon subspace is given by the expressionlij 

'P = Y.9i'j\^'^><^j\- (9) 

The expansion coefficients for the residues in real space 
as a function of the vector distance are given in Table ||. 
The coefficients for the transverse structure factor are 
given in Ref. 5. 
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We can now estimate the multi-magnon contribution to 
neutron scattering by simply subtracting A{q) from the 
equal-time correlation S{q). To get to the Heisenberg 
limit a series extrapolation is needed. Since the series for 
the multi-magnon weights are reduced by two terms (the 
first two being zero) compared to the parent series, one 
might suspect it would be better to extrapolate the series 
for the total cross section and the single-magnon contri- 
bution and take differences; we have used both methods 
to estimate the multi-magnon weights, carrying out the 
extrapolations by direct Pade approximants. In Fig. |^ 
results are presented for both the multi-magnon spectral 
weight as well as the ratio of the multi-magnon weight 
to the total spectral weight, along several lines in the 
Brillouin zone. We see that the multi-magnon contri- 
bution is particularly large near (tt, 0) (where, unfor- 
tunately, the extrapolation uncertainties are largest as 
well) and amounts to roughly a quarter of the total spec- 
tral weight. The spin- wave calculations of Igarashi and 
Watab(£3 yield roughly twice as much spectral weight 
in the multi-magnon excitations as we find by expan- 
sions about the Heisenhcrg-Ising model; however, given 
Igarashi's later remarkgj about the incorrect treatment 
of umklapp processes in that work, we do not view the 
discrepancy as a reason to question the accuracy of our 
results. |— I 

Recently, StringariliJ has developed general bounds 
and sum rules for single-magnon and multi-magnon spec- 
tra at special wavevectors. Near q = (0, 0) the sin- 
gle magnon spectral weight vanishes linearly while the 
multi-magnon weight vanishes quadratically; and near 
q = (tt, 7r), the single-magnon spectral weight diverges 
as Ktt, tt) — q|~^ while the multi-magnon weight goes to 
a constant. Our results are consistent with all of these 
requirements. 

In summary, series expansions and extrapolations have 
been carried out for dynamic properties of the S — 1/2 
square lattice Heisenberg- Ising model. Our numerical re- 
sults indicate that the Heisenberg model spin-wave spec- 
trum is close to but noticeably different from a uniformly 
renormalized classical spectrum. In addition, the single- 
magnon and multi-magnon spectral weights have been 
estimated throughout the Brillouin zone. In light of our 
numerical results, it would be interesting to examine the 
inelastic neutron scattering data on the antiferromag- 
netic parent compounds of the cuprate superconductors 
to look for the multi-magnon excitations. 
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FIG. 1. Spin-wave spectrum for the Heisenberg model, in 
units of J, along three lines in reciprocal space. The solid cir- 
cles with error bars are the results from the series expansions; 
the solid line is the classical spin-wave spectrum multiplied 
by an overall factor Zc ~ 1.18. 

FIG. 2. Heisenberg model multi-magnon spectral weight 
(crosses, narrowest error bars) and the ratio of the 
multi-magnon spectral weight to the total spectral weight as 
determined by extrapolations for the multi-magnon weight 
series (solid squares, widest error bars) and by the difference 
of extrapolations for the total weight and the single-magnon 
weight (diamonds, intermediate width error bars). 

TABLE I. Effective Hamiltonian for the Heisenberg-Ising 
model elementary excitations in real space, up to an overall 
factor of four. The dispersion in reciprocal space is found 
by summing all of the given real-space series with a factor 
(l/4)[cos(gir-i, + qyVy) + cos{qxra: - QyVy) + cos(q^r-y + qyr^) 



+ cos{qxry 


— Qyrx)], and then dividing by four. 


r 


Series 


(0,0) 


8 - 0.666662a^ + 0.6643520"* 




-0.292737a'' + 0.201076q* - 0.1774460^° 


(1,1) 


-4a^ + 1.222222Q'' - 0.5417560** 




-H0.513359a* - 0.430609q^° 


(2,0) 


-2a^ + O.IIHIIq" - 0.290209q'' 




-H0.255912Q* - 0.331992q1° 


(2,2) 


-0.2916670^ - 0.071979q'' + 0.060564a^ 




-0.133137a^° 


(3,1) 


-0.3888890^ - 0.177120q^ + 0.143711q* 




-0.220432q^° 


(3,3) 


-0.072627a*' + 0.002003q** - 0.029703q^° 
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(4.0) -0.0486110" - 0.074359a® + 0.021335a' 

-0.055154a"' 

(4.2) -0.108941a® - 0.013755a* - 0.050071a 

(4.4) -0.016444a* - 0.010852a^° 

(5.1) -0.043576a® - 0.030875a* - 0.034415a 

(5.3) -0.026311a* - 0.021359a^° 

(5.5) -0.005482a"' 

(6.0) -0.003631a® - 0.009438a* - 0.011028a 

(6.2) -0.013155a* - 0.017717a^° 

(6.4) -0.009137a^° 

(7.1) -0.003759a* - 0.009967a^° 

(7.3) -0.005221a" 

(8.0) -0.000235a* - 0.001596a" 

(8.2) -0.001958a" 

(9.1) -0.000435a" 
(10,0) -0.000022a" 



TABLE II. Single-magnon spectral weight series in real 
space. To evaluate the residue A(q) carry out the sum de- 
scribed in the preceeding table caption. 



r 


Series 


0,0) 


0.5 - 0.041667a^ + 0.011685a* 




-0.030642a® + 0.024677a* 


1,0) 


-0.666667a + 0.110185a^ - 0.130751a^ 




-|-0.117059a^ - 0.153300a^ 




0.194444a2 + 0.084876a" - 0.035514a® 




+0.092202a* 


2,0) 


0.097222a2 + 0.072647a" - 0.015601a® 




+0.057647a* 


2,1) 


-0.216667a^ - 0.096930a^ + 0.055191a'' 




-0.15U90a^ 


2,2) 


0.065365a" + 0.008393a® + 0.030857a* 


3,0) 


-0.036111o-^ - 0.049374a^ + 0.006892a'' 




-0.045770a^ 


3,1) 


0.087153a" + 0.034750a® + 0.054776a* 


3,2) 


-0.074005a^ - 0.018980a'' - 0.049672a^ 


3,3) 


0.018793a® + 0.015430a* 


4,0) 


0.010894a" + 0.018817a® + 0.018378a* 


4,1) 


-0.037003a^ - 0.028024a'' - 0.036609a^ 


4,2) 


0.028190a® + 0.029517a* 


4,3) 


-0.023254a'' - 0.021256a^ 


4,4) 


0.007614a* 


5,0) 


-0.003700a^ - 0.009765a'' - 0.011880a^ 


5,1) 


0.011276a® + 0.021890a* 


5,2) 


-0.013953a'' - 0.018825a'' 


5,3) 


0.012182a* 


5,4) 


-0.009425a^ 


6,0) 


0.000940a® + 0.004772a* 


6,1) 


-0.004651a'' - 0.011938a^ 


6,2) 


0.006091a* 


6,3) 


-0.006284a^ 


7,0) 


-0.000332a'' - 0.002181a'' 


7,1) 


0.001740a* 


7,2) 


-0.002693a^ 


8,0) 


0.000109a* 


8,1) 


-0.000673a" 


9,0) 


-0.000037a" 
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